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Abstract 

This thesis presents an algorithm suitable for numerical analysis of cryogenic refrigeration 
systems. The need for such an algorithm arises from the absence of effective analysis tools in 
the cryogenic field. 

Typically, preliminary design of a cryogenic system commences with a number of decoupling 
assumptions with regard to the process variables of heat and work transfer (e.g. work input rate, 
heat loading rates) and state variables (pinch points, momentum losses). These assumptions are 
made to facilitate direct and simplified solution calculations. However, making preliminary 
performance estimations minimizes the effect of component interactions which is inconsistent 
with the intent of "analysis". A more useful design and analysis tool is one in which no restrictions 
are applied to the system - interactions become absolutely coupled and governed by the equi- 
librium state variables. Such a model would require consideration of hardware specifications 
and performance data and information with respect to the thermal environment. Model output 
would consist of the independent thermodynamic state variables from which process variables 
and performance parameters may be computed. This model will have a framework compatible 
for numerical solution on a digital computer so that it may be interfaced with graphic symbology 
for user interaction. 

This algorithm approaches cryogenic problems in a highly-coupled state -dependent manner. 
The framework for this algorithm revolves around the revolutionary thermodynamic solution 
technique for Computer Aided Thermodynamics (CAT) developed by Dr. Gilberto Russo at 
the Massachusetts Institute of Technology. Fundamental differences exist between the Control 
Volume (CV) algorithm and CAT, which will be discussed where appropriate. 

This thesis presents the algorithm in detail with respect to the thermodynamic principles, 
state-variable management, device models, system networking, solution structure, numerical 
reduction, model flexibility and generality and potential pitfalls. 

Thesis Supervisor: Professor Joseph L. Smith, JR. 

Title: Professor of Mechanical Engineering 

Thesis Supervisor: Dr. Gilberto Russo 

Title: Lecturer, Department of Mechanical Engineering 
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1 Introduction 



1.1 Background 

Much design work has been done in the cryogenic engineering field and, accord- 
ingly, there are numerous design methods that have proven useful. Design, in an engi- 
neering context, means to model physical behavior and assign value to problem variables 
in such a manner as to achieve a desired engineering effect. A "good" design 
distinguishes itself from a "bad" design by whether or not the variables and parameters 
were assigned value in an optimal sense. Optimal is relative, but usually implies maxi- 
mum desired effect at minimal resource cost. Efficiencies and dimensional analysis vari- 
ables are typical ways in which engineers evaluate the optimization of a particular design. 
Thus, good design requires quality optimization. Quality optimization requires a thorough 
understanding of the problem variables, or, more succinctly, good analysis. 

Analysis means understanding the interactive behavior of the relevant physical vari- 
ables. This is the reason good design work is supported by a sound foundation of analyti- 
cal study. Good analysis tools are not in abundance in the cryogenic field, which forms 
the basis for this study: to develop an analysis algorithm suitable for application to 
cryogenic engineering systems. A typical analysis method for a closed loop thermody- 
namic system starts with a definition of the physical system in terms of models for 
devices (turbines, boilers, etc.) and for their interaction with the environment (Q^,, W out , 
T^etc.) . To determine the system performance, a number of assumptions must be made 
to reduce the vast number of unknowns. Such assumptions include the thermodynamic 
state at the at the inlet/outlet of a device (e.g. saturated, superheated, quality), a ratio of 
relevant pressures (AP/P, P H /P L ), or a temperature change ratio (T H /T L , pinch points). The 
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effect of making these assumptions is to decouple or minimize the interaction of state 
variables that govern system behavior; the behavioral patterns of the system become 
locked out. 

A robust analysis is one which minimizes global assumptions allowing the process 
variables to result directly from the highly coupled interaction of equilibrium state vari- 
ables. The analysis method could be likened to a black box which takes inputs of compo- 
nent specifications and yields as output the value of the system’s equilibrium state 
variables. It is then a simple extension to transform the state variables into more 
significant performance parameters. Thus system modelling is achieved by allowing the 
system to operate in its unrestricted condition compatible with the hardware imposed. 
Optimization (and hence design) becomes a matter of selecting hardware that achieves 
the best overall system performance, as defined by the designer. The method which is the 
heart of this paper pursues this strategy. 

Analysis is difficult. It requires simultaneous segregation and synthesis of the prob- 
lem variables to establish process characteristics or system behavior. Many engineers 
incorrectly use the word analysis when referring to design calculations. In fact, analysis 
in its crudest (yet often practiced) form is executed by conducting a series of design cal- 
culations with differing input values and observing the output or performance values. 

This is analogous to the experimenter plotting the results of a run without first calculating 
the analytical prediction. Such an analysis method does the segregation while loosing 
sight of the synthesis. This type of investigation is more appropriately termed design 
analysis. 
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1.2 Basis for Method 



In this work, the methodology for analysis of a thermodynamic system consists of 
idealization of the physical problem, generation of the mathematics modelling the prob- 
lem, solution of the mathematical representation and interpretation of the results. The 
aggregate of these tasks is rather complex. The method presented here reduces the 
complexity by using simple concepts to construct a complex representation of thermody- 
namic behavior that will be managed by a computer. This is essentially what is done in a 
finite element analysis method: building a complex representation to be managed by a 
computer from almost trivially simple engineering concepts. 

The great utility of a digital computer is that it can be programmed to perform the 
equation(s) generation and mathematical solution (through an intelligent assembly of the 
"simple parts"), leaving the idealization and interpretation to the user (although some of 
today’s modem routines will assist in the interpretation process by a convenient assembly 
of the output data). The method to be developed will allow a graphic, symbolic input of 
the problem, thus further assisting the user in the analysis of alternatives. The ability of a 
computer to manage large quantities of data (specifically, unknown variables) eliminates 
the temptation to make decoupling assumptions about problem variables that reduces 
their number to an (algebraically) manageable level. It is this special feature that makes 
this method unique: elimination of decoupling assumptions with respect to the problem 
variables causes the variables to become mathematically unconstrained and therefore 
physically interactive through the interconnective matching conditions. 

In applying a computer to perform thermodynamic analysis, the method employed 
should be general enough to allow analysis of a wide variety of complex multiple-element 
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thermodynamic systems. In such an employment, it should use simple models (vis. equa- 
tion sets) to represent the composite system and the same solution procedure to solve for 
the requisite data. This was the inspiration behind Computer Aided Thermodynamics 
(CAT) as developed by Dr. Gilberto Russo at the Massachusetts Institute of Technology. 
CAT is structurally similar to a finite element analysis except that it is applied to a collec- 
tion of discrete elements (that may undergo different types of interactions) rather than a 
homogeneous, continuous system (undergoing a single interaction). The great success of 
CAT has served as a backdrop for the Control Volume (CV) algorithm: adapting thermo- 
dynamic analysis to take advantage of computer power by formulating a new method of 
analysis. 
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2 Foundations of Control Volume Method 

2.1 Overview 

The control volume analysis method is developed for implementation on a digital 
computer. The basis for the method is modelling of complex systems using simple build- 
ing elements. These elements come in two broad classifications: storage elements and 
interconnective elements. Each element has associated with it a set of constitutive 
relations that describe the element’s behavior and are defined in a problem independent 
manner. 

The complete physical system is represented by all the constitutive relations of all 
the elements making up the system specified for the particular arrangement of elements. 
Once the elemental constitutive relations become specified for the system (through a 
change of base or variable transformation), they become the basis for the residual rela- 
tions of the system. Residual relations are used to define the steady state behavior of the 
entire system. The specificity comes from invocation of the system topology. The system 
topology is used to perform the change of base that generates the problem specificity 
from the "n" sets of constitutive relations of the "n" elements of the mesh. This accom- 
plishes the ’modelling and equation generation’ steps in this computer-aided thermody- 
namic analysis. Solution is found by using the residual functions in a global equilibrium 
equation: 



[rooi {A*} = {/?(*,)} CD 

[K‘(x)] is a matrix of first derivatives of the residual relations {R(Xj)}with respect to the 
independent variables {xjof the system (i.e. the Jacobian matrix). Successive modifica- 
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tions are made by imcrementing the current value of {x} by (Ax j which will drive 
[K‘(Xj)] (Ax) to zero (or less than some appropriate e). [K'txJKAxJ^O satisfies the 
constitutive relations of the system. The independent variable vector { x ( } for which 
(R(x)} goes numerically close to zero is the equilibrium state of the system and no fur- 
ther increment in x is required. 

[K‘(Xi)] is defined as the System Tangent Stiffness Matrix 1 . When [K‘(Xj)] is 
assigned a numerical value using (Xj), [K‘(Xj)] is the linear approximation to the system’s 
residual functions (R(Xj)}. This linearization is necessary due to the non-linear nature of 
some of the system’s equations. The System Tangent Stiffness Matrix represents an (n-1) 
dimensional tangent surface to a (n) dimensional function, hence the tag "tangent". The 
modifier stiffness comes from the congruence between equation (1) and the equilibrium 
equation for a structural finite element analysis system where the constitutive matrix is 
often referred to as the "stiffness" matrix. 

This paper describes the physical idealization of the thermodynamic problem, gen- 
eration of the mathematics describing the problem and solution procedure. Interpretation 
of results is not considered here, since that is a matter of recasting equilibrium state 
variables into more pragmatic (or new) performance representations (COP, rj^, etc.) that 
may be defined in the post-processing routines. 



1 "A New Methodology of Computer-Aided Thermodynamics", Gilberto Russo [doctoral 
thesis], Department of Nuclear Engineering, Massachusetts Institute of Technology, January, 



1987. 
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2.2 Compatibility 



The algebraic uniqueness and simplicity of each storage element lead one to believe 
that system simulation is now just a matter of connecting together elements of interest 
and performing the system integration and reduction. However, the issue of compatibility 
must be addressed with respect to the interaction of the engineering devices (which are 
independently defined). In other words, does the ’whole equal the sum of the pans’. 

This problem commonly arises during the design of thermal power systems with 
regard to "matching" of components. For example, does the mass-flow/pressure ratio 
characteristic of a compressor match that of the turbine it is to operate with? If not, it may 
be impossible to achieve a steady state operation for the system in the absence of spe- 
cially designed controls. Such a system (without controls) may oscillate (its performance 
and state variables) about a stable equilibrium (i.e. hunting), one component may drive 
the other to operate at an inefficient operating point. The stability of the operating state 
will require further consideration of the dynamics of the system. 

Accordingly, this methodology investigates the steady state operation of systems for 
which performance parameters (constants or functions) are assigned. Furthermore, the 
relations and models used to describe component behavior must be consistent with physi- 
cal limitations of a system. For example, it is incompatible to specify pressure rises for 
both a compressor and an expander which operate in the same fluid circuit since this 
would be to over-constrain the problem. Only one pressure ratio can be specified while 
the other pressure ratio results as a consequence of the interaction of the pressure and 
mass flow variables in the system. This compatibility is necessary to ensure a steady state 
can be attained. 
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To accommodate the compatibility issue during the initial stage of the development 
of the method, pressure ratios will be defined only for compression devices. All expan- 
sion devices (turbines, throttles, etc.) will be void of any pressure characteristic; the equi- 
librium pressures will be established by the interaction of the system elements. 

2.3 Physical Idealization 

The physical system is modelled using its constitutive elements . Figure 1 is a repre- 
sentation of a closed loop thermodynamic circuit operating as a refrigeration cycle. While 
heat and work interactions are an intimate part of the complete operating system, these 
external interactions are determined after the equilibrium behavior of the system is estab- 
lished. The work or heat interaction for any element is readily determined by simply 
knowing the inlet and outlet state. This is infact the approach used here: work or heat 
transfers are determined after the value of the equilibrium state variables are known by 
writing the 1st Law of thermodynamics in the post processing routines. 

The system is the collection of control volumes that contains all components of 
interest and the environment with which the components interact. A valid thermal analy- 
sis of the component interactions can only be done by analyzing a uniquely independent 
set of control volumes. This system has four distinct components 2 . Accordingly, a set of 
four independent control volumes must be developed. A selection of any three of the four 
control volumes shown in Figure 1 plus a control volume encompassing the entire system 
(including the environment, which is part of the system) is a valid selection. However, 



2 For simplicity in this development, interconnective piping and flow passages are assumed 
to be lossless. Lossy flow in interconnecting piping may be incorporated into the framework 
upon extension of the applicable constitutive relations. 
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the control volume methodology is based on using simple construction elements to build 
the total system model, therefore the four component control volumes (including their 
external interactions) idealize the system. 

In choosing the relations that will be sufficient to describe the system being ideal- 
ized, careful thought must be employed. If an engineering system is not completely speci- 
fied in terms of the relevant behavioral relations, then the number of additional arbitrary 
parameters that must be specified will equal the number of degrees of uncertainty left in 
the mathematical representation of the physical system. On the other hand, if the system 
is overspecified, ambiguous and contradictory results will precipitate (i.e. over- 
constrained, physical incompatibility). There must be a compatibility between the physi- 
cal and mathematical problem: the number of independent unknowns must equal the 
number of independent relations. One of the obstacles of this effort was in finding this 
delicate balance between the mathematical and physical requirements. 
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Figure 1. : Conventional representation for analysis of a thermal circuit 



2.3.1 Local and Global Variables 

Local variables are independently defined within each storage element and are 
used in the constitutive relations of each storage element. When elements are 
assembled together to form a system, the local variables are now transformed into a 
variable system defined for the complete physical system. This system of variables is 
called the global variables system and it is defined only for a specific system based on 
the specific interconnections between storage elements. 

2.3.2 Storage Elements 

A storage element is in a pragmatic sense a self-contained standard engineering 
component such as a turbine, heat exchanger or valve that is defined on an elemental 
basis. Storage elements transfer energy, entropy, and mass to achieve some desired 
engineering effect. Accordingly, the 1st Law of Thermodynamics, continuity require- 
ments, and an efficiency (irreversibility parameter) relation are chosen to describe the 
constitutive behavior of the element. 

For devices which exchange work with their environment, the classical defini- 
tions of component isentropic efficiency serves as the irreversibilty parameter. Such a 
definition reflects 2nd Law irreversibilities while incorporating component specific 
"test-stand" data. 

The irreversibility relation for devices which exchange heat (either internally or 
with the environment) is from the heat transfer rate equation. The irreversibility infor- 
mation is contained in the overall heat transfer coefficient (U) and the available heat 
transfer area (A) which are used in the rate equation to relate fluid stream 
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temperatures. A second irreversibility parameter is from the fluid friction losses. Thus, 
relations defining the pressure drop through a heat exchanger together with the heat 
transfer rate equation define the irreversibility (entropy generation). This model may 
seem to make the methodology more complex, but it achieves the goal of incorporat- 
ing hardware performance data into the elemental models. 

The equations used to specify the irreversibilties are called the characteristic 
relations since they reflect component’s characteristic behavior. 

The working fluid constitutes an additional constitutive relation for each fluid 
state. The type of working fluid employed should be arbitrary in a useful analysis 
method. The physical properties of the fluid should be expressible as a function of 
independent variables. This has been accomplished for the actual fluid states, but for 
the isentropic outlet states (work transfer devices), working fluid properties are 
expressed as a function of the independent pressure variable and the value of the isen- 
trope. That is, h out =h(P out , s J. 

If a system uses a single working fluid, then all states are calculated using the 
same property relation. On the other hand, if circuit contains multiple, independent 
fluid loops each employing different fluids(e.g. a nitrogen pre-cooler for a liquid 
helium plant), then a fluid designation must be associated with each loop. In other 
words, each fluid loop would have a working fluid constitutive relation applicable to 
that fluid. This is more efficient than defining the fluid constitutive behavior for every 
element. 
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In summary, each storage element has constitutive relations grouped categori- 
cally by 1st Law, characteristic relations, mass continuity and working fluid constitu- 
tive behavior. 

2.3.3 Interconnective Elements 

The interconnective elements model (through their constitutive relations) the 
linking of various storage elements which constitute the system. As individual compo- 
nents are assembled with their associated elemental constitutive relations and local 
variables, a transformation is necessary to shift from the local variables system to the 
global variables system. Such a transformation is accomplished by the interconnective 
element’s constitutive relations. By transforming from a local to a global variables 
system, topological information is implicitly contained in the new variables system. 
The interconnectivity requirements (i. e. the constitutive relations of the interconnec- 
tions) enforce on the system the equality of independent local variables at a nodal con- 
nection between storage elements. Since these local variables are equal, they assume 
the global variable label at the node. This ensures component interaction. For 
example, in Figure 1, the outlet pressure from the compressor must equal the inlet 
pressure to the high temperature heat exchanger. The two local variables are equal and 
are assigned the global pressure variable. Similar operations are performed on all the 
independent variables (Chapter 3). When the transformation is complete, the number 
of global variables equals the number of local variables less the number of intercon- 
nection constitutive relations. 

2.4 Equation Generation 

The constitutive relations of the system are expressed as the { R(x s ) } vector which is 
derived from the constitutive relations of the elements and the topology of the system. 
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Equation generation refers in this context to the generation of the [K‘(x*)] matrix and the 
{ R(x ( ) ) vector. These quantities are generated through recursive relations which draw on 
topological information to generate the appropriate elements of both [K‘(x)j and {R(x)}. 
The methods for generating these relations are discussed in further detail in Chapter 6. 

2.5 Solution 

The solution technique, briefed in Section 2.1 is a Newton’s method tailored for this 
application and will be reviewed in further detail in Chapter 7. 
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3 Variables and Parameters Management 
3.1 Selection of Independent Variables 

Section 2.2 made indirect reference to the properties of temperature and pressure as 
the independent thermodynamic variables of this method. In fact, a significant portion of 
the CV algorithm development was performed with temperature and pressure as the inde- 
pendent variables. However, temperature and pressure are intensive thermodynamic prop- 
erties. Two intensive properties can fix the thermodynamic state in a homogeneous, 
single phase system only. In generalizing this algorithm to an operating refrigeration 
cycle, determination of properties and states in a binary or two phase system is necessary, 
since much of the refrigeration effect results from the latent heat. Therefore, to lend gen- 
erality to the algorithm, pressure and specific volume were selected as two independent 
variables. 

For shaft work machinery(turbines, compressors, etc), three equilibrium states are 
of interest. The equilibrium inlet state, the isentropic outlet state and the actual outlet 
state. The isentropic oudet state can be determined as a function of the inlet entropy value 
and outlet pressure value: h 2s =h(P 2 ,s 1 ), where s^s^P^v,). 

In addition to pressure and specific volume, mass flow rate has also been selected as 
an independent variable. Mass flux scales the system’s energy and entropy interactions. 
Mass flow rate is an essential variable when devices such as splitting valves and mix- 
ing/recombination valves are employed. 
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These three variables (pressure, specific volume, mass flow rate) taken together 
with the interconnectivity relations (between devices) ensure full thermodynamic match- 
ing at the "nodes" between storage elements. Thus, a complete engineering link is made 
between two devices exchanging fluid energy. 

3.2 Dependent Variables 

The dependent state variables are all other thermodynamic properties: 

h(P,v) 

h s (P,s(P,v)) 

s(P,v) 

T(P,v) 

Work and heat transfer are also system unknowns, although not variables perse. They are 
determined by the equilibrium values of the global variables. The 1st Law residual func- 
tions and [K l ] entries for elements which exchange heat or work with the environment are 
not incorporated as part of equation (1) framework. Such relations are assigned to a 
post-processing routine. In doing this, no specificity of the [K‘(x)]{Ax}={R(x)} equilib- 
rium equation is lost since elimination of a 1st Law relation also eliminates either a work 
transfer or heat transfer unknown. 

Determination of the equilibrium state variables { x ) completely defines the state of 
the system, from which heat and work interactions can be determined. 
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3.3 System Parameters 

System parameters are user defined as part of the algorithm initialization where the 
specifications of component data are imported for use in the applicable constitutive rela- 
tions. Such parameters include component efficiencies Cn t ,T|c)> friction factors (f), heat 
transfer coefficients (U), heat transfer areas (A), and expansion/compression ratios (r E , r c ). 
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4 Device Models for Control Volume Analysis 

4.1 General 

A set of constitutive relations for the individual control volumes has been devel- 
oped which fit into a standardized framework based on the engineering principles dis- 
cussed in Chapter 2. The constitutive relations are grouped categorically into 1st Law, 
continuity, interconnectivity and characteristic relationships (characteristic relations 
describe the non-ideal behavior of the selected control volume in view of the Second Law 
irreversibilities that degrade the system performance from the ideal case). The categorical 
standardization of the relations for different types of engineering devices is necessary to 
facilitate the assembly of the elements of [K*( { x } )] { Ax } = { R( { x } ) } via numerical meth- 
ods. 

4.2 Engineering Devices 

The following devices have been selected for implementation in the control volume 
analysis: 

Regenerative Heat Exchangers 

Ambient Heat Exchangers 

Compressors 

Expanders 

Expansion Valves 

While this list is not all-inclusive, it is representative of typical devices found in 
cryogenic refrigeration systems and thermal power systems in general. More complicated 
devices such as Stirling engines are not covered here but may be included upon extension 
of the constitutive framework discussed in Chapter 2. 
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4.3 Component Constitutive Relations 

Symbology for the expander is shown in Figure 2. 



1 st Law: 



^A(P l ,v l )-m 2 h 2 (P 2 ,v 2 ) = W E 



Characteristic relations: 



(/*i /t 2 ) h is ) — o 

h i-^ = 0 

The first equation is the definition for isentropic efficiency. The second equation is a con- 
sequence of the recursion used to generate the efficiency equation: h2s can not be numer- 
ically generated without generating his. This is because the recursive relations use 
topological information to generate the appropriate functions. Topological information 
can discriminate between inlet and outlet states but not between isentropic and actual 
states. The recursion generates an isentropic inlet state which is identical to the actual 
inlet state; the second equation merely sets the two states equal. The third equation 
equates the inlet entropy to the outlet entropy for an isentropic process( which must be 
known to use the efficiency definition). is used in a fluid constitutive relation to deter- 
mine h 2s =h(P 2 , s*). 

continuity: m,-m 2 = 0 

working fluid: h^-h (P 2 , %) = 0 
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Fig 2. : Expansion Engine Schematic 




Fig 3. : Compression Engine Schematic 
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The working fluid constitutive relation defines the outlet enthalpy for the isentropic pro- 



cess. 



The nine unknowns in these equations (two mass fluxes, two pressures, two specific 
volumes, an isentropic enthalpy and entropy and work rate) are represented by only six 
equations. Since each control volume is a basic building element of the CV algorithm, 
each individual control volume should be algebraically independent. Therefore three 
additional equations are necessary - these equations come from the interconnectivity 
requirements for pressure, specific volume and mass rate. Each (2-port) control volume 
will have associated with it three interconnectivity relations (from an adjacent intercon- 
nective element) which add to the other constitutive relations to ensure algebraic unique- 
ness. 



The applicable constitutive relations can also be developed for the compressor 
engine, which is schematically shown in Figure 3. 

1st Law: 

m A(P l ,v l )-m 2 h 2 (P 2 ,v 2 ) = W c 



Characteristic relations: 

thermal characteristic: 



T\c( h \-h?)-(K-K) = Q 
h x ~! »,, = () 
s*- Si =0 
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pressure characteristic: 



p i~ r c p 2 = 0 



mass continuity: 



m l — m 2 = 0 



working fluid: 



^ 2 j ^ (P $2s) ~ 0 



These equations are similar to those for the expander, accept that a pressure rise relation 
(which sets the system pressure levels) has been incorporated. 

A schematic for the regenerative heat exchanger is shown in Figure 4. 

1st Law: 



m x h x (P v v x ) -th 2 h 2 (P 2 , v 2 ) + m 2 h 2 (p v v s ) - m A(p 4 , v 4 ) = 0 



thermal characteristic: 



m x h x — m 2 h 2 — U R A R lnj 



(T X -T A )-{T 2 -T 2 ) \ 
(T X -T 4 )/(T 2 -T 2 ) J 



= 0 



pressure characteristic: 



2 " 






m 2 i(v l + v 2 ) = 0 
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Fig 4. : Regenerative Heat Exchanger Schematic 





DISCONTINUOUS AT a = AT b 



Fig 5. : Temperature profiles for which the rate equation is invalid 



p,-p*~ 



16/7 ) , 

^J / 3(V3 + V4) = 0 

continuity: 

m { -m 2 = 0 m 3 -m 4 = 0 

The temperature characteristic is a form of the rate equation. Stream 1-2 is assumed to be 
the warm stream and stream 3-4 is the cold stream. It is a valid representation of the tem- 
perature profile through a heat exchanger for profiles that exhibit a continuous, exponen- 
tial temperature difference along the flow path. Figure 5 shows temperature profiles for 
which the rate equation is not valid: discontinuity in the profile or invariant profile. In 
instances where a phase change occurs, two separate rate equations should be used, one 
for each side of the discontinuity. For an invariant profile, the rate equation is valid by 
taking a limit when ATa/ATb goes to 0/0. "UA" is the component heat transfer coefficient 
and heat transfer area. 

Two pressure relations are given here for the friction loss phenomena. For model 
development, the pressure drop is simply represented by a flow duct with characteristic 
hydraulic diameter, flow path and friction loss factor, which should be known apriori for 
a given component. This is a very generic representation of pressure drop due to friction. 
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The segregation of the characteristic equations for this device into temperature and 
pressure relations does not imply the independence of heat and momentum transfer. On 
the contrary, they are intimately coupled, which would be reflected by their respective 
values for "f ’ and "UA". 

An ambient heat exchanger is shown schematically in Figure 6. The constitutive 
relations for an ambient heat exchanger are very similar to those for the regenerative heat 
exchanger with the second flow stream being replaced by a quantity of heat being 
rejected to/ received from the environment. 

1st Law: 



m x h x {P v v,) -m 2 h 2 (P 2 , v 2 ) = Q 



temperature characteristic: 




pressure characteristic: 




continuity: 



m x -m 2 = 0 



33 






i 




Figure 6. : Ambient Heat Exchanger Schematic 




Figure 7. : Expansion (Throttling) Valve 
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The relations for the ambient heat exchanger are of the same form as those for the regen- 
erative heat exchanger. A subtle difference exists in that the ambient environment is rep- 
resented by a temperature reservoir T„. At steady state or for very large reservoirs ( i.e. 
^reservoir ^sys.toi)’ is invariant. T^yand U^A^ are also defined during initialization. 



The last device which is commonly found in cryogenic systems is the throttling or 
expander valve, shown schematically in Figure 7. 

1st Law: 



m l h l (P v v l )-m 2 h 2 (P 2 ,v 2 ) = 0 



characteristic relation: 



h\ ~ h 2 = 0 

continuity: 

m i -m 2 = 0 

The principle characteristic of the expansion valve is that the entering and exiting enthal- 
pies are equal (although the enthalpy through the valve is not constant). Note that the 
pressure drop across the valve is not specified as a characteristic relation. In order for the 
existence of this valve to make physical sense in a refrigeration system, it must operate in 
conjunction with a pressure elevating device (compressor). Placing a required pressure 
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ratio on the expansion valve would place a redundant overconstraint on the system, since 
the compressor already sets the system’s principal operating pressure levels (see Section 
2.2:Compatibility). 

4.4 Interconnectivity Relations 

The interconnectivity relations are mathematically trivial but form the critical link 
to the algorithm. By equilibrating the independent thermodynamic variables between the 
outlet of one device and the inlet of a connected device, thermodynamic matching is 
imposed on the system. Now, all elements behave in a highly-coupled interactive manner. 

Interconnectivity places three relations at each component interface: pressure, spe- 
cific volume and mass flux. Equilibration of pressure and specific volume ensure a ther- 
modynamic match of the working fluid and equilibration of mass flux closes the 
continuity loop around a closed system. Typically, interconnectivity relations will have 
form 

P 2 -P 2 = 0 v 2 -v 3 m 2 -m 2 = 0 

4.5 Summary 

The device models presented here are self-contained, problem-independent equation 
sets. The next task is to shape the constitutive relations so that they may be compatible for 
assembly of the [K r ] matrix and the {R} vector. 
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5 Topology 



5.1 Introduction 

Prior to discussing the assembly of the [K‘] matrix and {R} vector, it is useful to 
describe how the topology is developed, since [K'J and {R} are assembled using the sys- 
tem’s topolgy. The topology establishes in mathematical terms which elements are con- 
nected and how they are connected in an engineering system. Different types of 
interactions yield different topological structures or networks. 



The most efficient numerical method of representing the system topology is with 
the topology or incidence matrix. Such a matrix is a tertiary state (1, 0 or -1) matrix that 
relates the connectivity of nodal variables. Topology matrices are commonly found in 
discrete network systems such as electrical network circuits exchanging current or struc- 
tural networks exchanging displacements. In the Control Volume methodology, topology 
is used to describe how elements exchange mass, energy and entropy. Hence, there is 
more than one interaction at each node. Figure 8 shows a typical structural network. The 
springs represent elements that store energy; the nodes represent equipotential points 
between storage elements. The topology matrix for this system is 



Top[i,j] = 



1 -1 

0 1 

0 0 

-1 0 

.-1 0 



0 0 

-1 0 

1 -1 

1 0 

0 1 j 
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Figure 8. : A structural network; kj represent elemental stiffnesses 
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The columns represent nodes (interconnective elements) and the rows represent storage 
elements. In the CV algorithm, topology is defined on two levels: thermal and mass 
exchange. Thermal topology describes elements in thermodynamic communication: stor- 
age elements exchanging heat, work, fluid energy and entropy. The mass exchange topol- 
ogy describes the arrangement of elements exchanging mass and momentum. 

Why is this done? It is principally a matter of mathematical convenience for the assembly 
of the appropriate residual functions and [K‘] functions. The thermal topology is used to 
assemble the 1 st Law and temperature characteristics, while the mass topology matrix is 
used to assemble the pressure characteristic and mass continuity relations. The distinction 
arises due to the presence of counterflow heat exchangers that have one 1st Law relation- 
ship and one temperature characteristic but two pressure loss functions and two mass flow 
relations. In such an arrangement, the two separate flow passages in a counterflow heat 
exchanger each constitute a separate pressure loss mechanism. But only one energy inter- 
action is present. 

The computer will generate these matrices by scanning screen. However, the ini- 
tialization routine will require two separate topological generations, one on the thermal 
interaction level and one on the mass transfer interaction level. 
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As an example of the application of the control volume toplogy, consider the 
Claude cycle shown in Figure 9. Five storage elements and six interconnective elements 
(nodes) are shown. This system will have five 1st Law relations and thermal interaction 
characteristics, but six momentum interactions since there are six distinct flow passages. 
The thermal topology matrix is 



[Top T ] = 



-1 0 

0 0 

1 - 1 

0 1 

.0 0 



0 0 

0 0 

0 1 

-1 0 

1 -1 



0 1 

1 -1 

-1 0 

0 0 

0 0 _ 



The mass and momentum exchange topology matrix is given here for the same 
Claude cycle, but the schematic difference is shown in Figure 10. The distinct flow cir- 
cuits of the regenerative heat exchanger are now modelled as separate entities such that 
the corresponding topology can be used to generate the pressure characteristic relations 
and the mass flux relations. 



[TopJ = 



-1 0 

0 0 

0 0 

1 -1 

0 1 

.0 0 



0 0 

0 0 

0 1 

0 0 

-1 0 

1 -1 



0 1 

1 -1 

-1 0 

0 0 

0 0 

0 0 _ 



The convention used in both topological representations is taken as flow into an element 
as positive and flow out of an element as negative. 
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Figure 9. : Thermal topology representation for a Claude cycle: 
six nodes, 5 storage elements. 
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Figure 10. : Mass topology representation for a Claude cycle: 
six nodes, six storage elements. 
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6 Generation of the Residual Relations 

6.1 Introduction 

Chapter’s 4 (elemental models) and 5 (topology) laid the mathematical foundation 
for the constitutive relations that serve as the basis for generating the [K'] matrix and the 
{ R) vector. The algebraic equation sets presented in Chapter 4 must be manipulated to 
conform to the numerical requirements of [K‘(x i )]{Ax}={R(x i )}. The [K‘] matrix is gener- 
ated from the elemental constitutive relations in two steps. Since [K‘] is a Jacobian of the 
system residual relations, the first derivatives (with respect to the global variables) of all 
the elemental constitutive relations must be computed and these first derivatives are the 
constitutive relations of the elements assigned to their proper [K l ] locations. This is done 
using the system topology information. 

6.2 Assembly of the (R(x)} Vector 

{ R } is the array used to determine the steady state operating conditions of the sys- 
tem. It is also assembled using the system topology information. 

The assembly of the {R} vector is much simpler than the assembly of the [K l ] 
matrix (Section 6.4). A subtle differences between the generation of { R} and [K‘] is that 
the residual relations (R(x)} are multiple function equations(e.g. a 1st Law has several 
m./i, arguments summed together). Hence, their recursive relations shall contain summa- 
tion operators to cycle through the nodal variables and use the topology to filter out 

entries that are not pertinent to the specific residual. On the other hand, the elements in 

aq 

the [Kt] matrix are usually single-argument entities (e.g.^-)-no summation over the range 
of nodal variables is made. 
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As an example, consider the 1st Law relation for the regenerative heat exchanger 
shown in Fig. 10: 



m x h x -m 2 h 2 + m 4 h 4 -m 5 h 5 = R 

The thermal topology array (row) for the element in question will contain (+1) for the 
flows (from nodes) entering the element, (-1) for flows (to nodes) exiting elements and 
(0) for nodes not connected to the element under consideration. The regenerator in Fig. 10 
corresponds to storage element 3, which corresponds to row 3 in the thermal topology 
matrix Top-r[3,j] shown on page 42. Thus, 

Top r [3,/j = [1 -10 1-1 0] 

R= I rhjhjTop T [3,j] 

j = i 

where 6 equals the number of nodes in the system. Expanding the summation through the 
appropriate values of ’j’ (and hence Top T [3,j]) yields 

m x h x -m 2 h 2 + m 4 h 4 -m s h 5 = R 

Similar combinatorial operations can be performed to generate all the residual functions. 
However, not all constitutive relations have such well regimented structures as the addi- 
tive 1st Law. In such cases, the numerical value of the topological arguments are used to 
modify the variable indexes and parameters where appropriate. This is apparent in the 
pressure characteristic constitutive relation (Section 4.3) where one pressure is multiplied 
by the pressure ratio and the other pressure is not. 
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Px-rf 2 = 0 



The topolgy value (1 or -1) can be used to ensure that the pressure ratio multiplier appears 
only with the outlet pressure. The details of the recursion relations used to develop the 
{ R) vector are left to Appendix 1. 

6.3 [K‘] Matrix Functions 

The analytical functions that form the entries in [K‘] matrix are the first derivatives 
(with respect to the global independent variables) of the residual functions: 



Taken literally, the residual functions R s must be known prior to computing Kf. This 
would be the case for analytical generation of [K‘], but the control volume method is 
numerically based for implementation on a computer. It would not be the best employ- 
ment of a computer’s potential to generate and store all the system’s residual relations 
and then compute their derivatives. It would be prudent to store the unspecified first 
derivative constitutive relations (which are standardized for each element type) and then 
form system specific constitutive relations from these standardized relations by combina- 
tion with the system topologies. This is precisely the what is done here. 

As discussed in Chapter 4, all the possible "R’s" are limited to 1st Law, thermal and 
pressure characteristics and mass continuity relations for a limited collection of elements 
(Sect. 4.1). Xj are the system global parameters, namely pressure, specific volume, mass 
and flow rate at each interconnection. As an example, a 1st Law relation for a regenera- 
tive heat exchanger is developed as follows: 
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R,(Xj) = I m h = 0 

j = i 

where p represents the number of fluid streams entering or leaving the control volume for 
device "i". The partials of R with respect to one set of system global parameters are: 

dR, dR : dR, 

dP/ dv/ drhj 



These derivative functions have the following standard form: 



a/?, 



p dh. 



a/?, 

a V/ 



p dh. 

S m - - ] 

y = i J dV; p 



dR t 

drhj 



I 



j = i 



h. 



where the right hand side of the equations contain only constitutive relations and global 
variables. The use of the topology matrices for the generation of the [K l ] terms is pres- 
ented in Appendix 1. 

As mentioned in Section 2.3.2, 1st Law relations with heat and work transfer are not 
incorporated into the K' framework, so no derivatives need be computed for W or Q 
terms. W And Q computation is assigned to post-processing. 

Finally a numerical (instead of an analytical) derivative is used to represent the line- 
arization of the rate equation (heat exchanger). 
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6.4 Assembly of the [K‘] Matrix 

6.4.1 Interaction Regions of the [K l ] Matrix 

The constitutive relations that form the [K‘] matrix described in Section 6.3 are 
the basis for the formulation of the [K 1 ] terms. 

To develop [K‘], it is necessary to visualize the matrix consisting of separate 
interaction regions each corresponding to a type of interaction as described in Section 

2.3.2 and Chapter 4: 1st Law, thermal characteristics, pressure characteristics, and 
mass continuity. 



[K‘(xj\ {x} = 



1st Law 

temp, characteristic 
Pres, characteristic 
. mass continuity . 



W 



This administrative subdivision of the [K‘] matrix is necessary due to the two differing 
types of interaction that lead to two distinct topology matrices: thermal and mass 
exchange. The vector {x} is the collection of the independent variables relevant to the 
problem: {P X ,P 2 , • ••,P n ,Vx,v 2 , ...,v„,m 1 ,m 2 , ...,m n }. The derivative functions pres- 
ented in Appendix 1 must now be uniquely assigned to their appropriate location in 
[K l ] using the topological information. 

Since there will only be three types of independent variables, each interaction 
region within the [K‘] matrix may be further subdivided (in an administrative sense) 
into smaller 2nd-order arrays. 
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(R i = 1st Law) 



[*'(*)] {x} = 



(. R , = therm, char.) 
(Ri = pres, char.) 
(Ri = mass cont.) 



dR, 




dRi 






- dv j- 


d Ri 




'dR' 






__dv jz 


dRi 




dR, 


ZPj- 




dvj_ 


dRi 




'dR,' 






.dvj. 



' dRt 
dm 
' dR, 

driij 

dRi 

dm 

'a*, 

dm , 



j-i 



]-> 



V n 

m. 



M. 



Each array denoted by [ ] represents the collection of 1st derivatives of a constitutive 
relation with respect to one of the independent variables. The six arrays corresponding 
to 1st Law and thermal characteristic are identical in size to the thermal topology 
matrix, used to assemble the arrays from the constitutive relations of the elements. 
Similarly, the six arrays corresponding to pressure characteristic and mass continuity 
use the mass topolgy matrix for function generation. All arrays have the same number 
of columns (since each column represents a node), which ensures the columns pro- 
perly align when the total [K‘] matrix is assembled. 

The justification for this arrangement within the [K 1 ] matrix is that it makes 
assembly of the matrix simpler. Once the respective toplogy matrices are established, 
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each storage element is queried for its applicable constitutive relations (which form 
{R}) and [K 1 ]. These relations are combined with the topology through recursive rela- 
tions to define, in this case, the appropriate [K l ] entry. 

6.4.2 Array Assignments 

Each array shall be labelled using the following convention for ease in identify- 
ing what interactions it represents, where it should be assigned in [K l ], and for esta- 
blishing the recursions used to produce the [K‘] entires. The first superscript "t" 
indicates (as before) that the functions to be developed consist of terms of the system 
tangent stiffness matrix. The second superscript indicates the interaction type to which 
this array refers(lst Law =1, therm. char.=T, pres. char-P, mass cont.=m. The first 
subscripts are the array location indices. The outer subscript indicates the respective 
variable of derivation ( P, v or m). As an example, consider 

Superscript "1" means that this array refers to a 1st Law relation. Subscript ”P" indi- 
cates that the functions inside are the partial derivatives with respect to the indepen- 
dent pressure variable. The indices [3,4] indicate that this particular entry corresponds 
to the third 1st Law relation and the fourth global pressure. This labelling has been 
developed so that each inner matrix can be developed independent of the next one 
using the recursive relations that are to be presented in the following paragraphs. 
Using this identification scheme. The [K l ] matrix has as arrays 
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[K [ ] = 



— 

L__J 

1 




m T P 






m p v n p * 




mi 



The panial derivative functions of the constitutive relations of Appendix 1 (dis- 
cussed in Section 6,3) are combined with the respective topology matrices to form the 
inner matrices, which form the global [K‘] matrix. Each row within an inner matrix 
corresponds to the derivative of a storage element’s constitutive relations. 

As an example, the recursive relations to generate the inner matrices for the 1 st 
law interactions are: 



, i dh , 

(*»), =m^].Top r [/. J l 
( K^m^loVrV.i] 

(AT^ = A,ToftM 



The recursive relations for all inner matrices are contained in Appendix 3. 

6.4.3 Summary 

The [K 1 ] matrix has been subdivided for assembly purposes into interaction 
regions and inner arrays corresponding to the global variable. Each inner array is gen- 
erated by a recursion relation (App. 2) standardized for each type of control volume. 
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Any inner matrix may be identically zero if the constitutive relations (for that array) 
are not functions of the variable used to compute the partial derivatives. Each [K 1 ] 
term is a function of constitutive relations of the elements. 

Thermal toplogy is used to generate 1st Law and thermal characteristic functions 
while the Mass topology is used to generate the pressure characteristic and mass conti- 
nuity functions. If a system is composed entirely of two-port elements, then the mass 
and thermal topology matrices will be identical. In the general case however, the 
index ranges for the two matrices will be different. 



51 



7 Chronology of CV Method 



7.1 Introduction 

The formulation and solution procedure has three steps: preprocessing, generation 
(system assembly) and reduction. 

7.2 Preprocessing 

The first step in the preprocessing operation is to establish the control volumes of 
interest in a closed loop system. During the initialization of each individual storage ele- 
ment local variables are assigned based on the type of element, i.e. a standard two-port 
element (expander) will have inlet variables P 1 , v 1 , m 1 and outlet variables 

P 2 ,v 2 ,m 2 . A four port element will have P l ,...,P 4 etc. The local variable assignment 

sets the stage for the local to global variable transformation (i.e. system specification) as 
discussed in Sections 2.3.1, 2.3.3 and 4.4. 

These transformations are conducted for all the independent variables (P and v as 
well) at the interconnecting nodes only. This local to global transform has two functions: 
(1) it reduces the number of independent state variables by a factor of two thereby making 
the system algebraically unique and (2) it impresses on the system the interactive require- 
ments of the individual elements. 

When the elements that constitute the system are selected and connected per desig- 
ner’s requirements, the interconnectivity is then known and the topology can be estab- 
lished. The topology is then used to generate the applicable matrices and vectors. 
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7.3 System Assembly 

The assembly revolves around the individual control volume as a basic construc- 
tion element. Once the algorithm determines the number of storage elements and connec- 
tive nodes (via the topology matrices), it will establish space (memory) sufficient to carry 
the requisite number of constitutive relations based on the characteristics of each control 
volume. The details of the assembly operation were discussed in Chapter 6. 

7.4 Reduction 

Reduction is the actual numerical solution of the [K t (x k )]{Ax k }={R(x k )}=0 system. 
It uses a simple Newton-Raphson method to successively iterate toward the equilibrium 
state, {x eq } (the index "k" represents the "kth" numerical value based on thermodynamic 
state {x k }). This method is a common numerical exercise and is discussed in detail in 
Appendix 4. 
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8 Conclusions 



This effort was the first step in an effort that is to develop new computer-based tech- 
niques for the simulation of cryogenic engineering plant performance. This work centered on 
cycle analysis of steady-flow cryogenic systems. Simulation is the extension of this analysis 
to more generalized operating scenarios. Analysis is indeed difficult, yet it is only the first 
step toward simulation. In analysis however, some important foundations can be established 
(or alternatively, bad foundations can be eliminated) which may be incorporated into simu- 
lation models. Such matters concern variable management, idealization of the physical prob- 
lem, mathematical formulation of the idealization and a solution algorithm. There are some 
fundamental obstacles which stand in the way of achieving a simulation model for a thermal 
power system. These obstacles arose during this effort conspicuously and inconspicuously- 
certain obstacles remain, some were overcome. The remarks that follow address in retrospect 
the extension of this work to future development and the problems in simulation in general. 

Most of the work that has been done in the simulation of thermal power systems has 
been done in the nuclear power plant field due to the critical need for understanding of plant 
behavior before embarking on construction. However, most simulation models are limited to 
a very select set of flow circuits. A methodology to perform analysis (or simulation) of arbi- 
trary thermal circuits has not been rigorously developed or disseminated. This paper pres- 
ented a method of analysis for arbitrary thermal power systems. 

Another fundamental obstacle in analysis and simulation of power systems lies in the 
historical practice of thermodynamics. The industrial age has seen the expansion of power 
sources (electrical, automotive) to all comers of the earth. Accordingly, much of the thermal 
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power engineering has been design oriented. As a result, thermodynamic models for pro- 
cesses have been tooled to yield design data based on some required performance specifica- 
tion: 1st Law results scale sizes of systems; efficiency definitions dictate what efficiency is 
required to achieve some net work; rate equations define area and heat transfer requirements 
for a specified heat transfer. Such models form the cornerstone of time-proven design aids. 
Design is the antithesis of analysis and simulation, and herein lies the difficulty for the analy- 
sis problem: a plethora of good design models, yet a dearth of useful analysis models. What 
is frequently done is the design models are "reverse-engineered": performance is determined 
based on specified component design data. The practical application of such an approach has 
shortcomings. 

Design models are basically lumped-parameter relations (whereas steady flow devices 
usually exhibit gradients. These lumped-parameter models are useful due to their simplicity 
which was one reason for implementing them in this method: use simple models to build 
complex systems. More complex elements will be handled in the same way by the processing 
and solution algorithm. However, as was seen in Chapter 6 (and in the appendices), even 
these simple models required significant algebraic and numerical manipulation to fit the into 
framework. Another shortcoming of the lumped-parameter formulation was that it restricted 
the class of problems to be studied (steady flow systems). To analyze off-design and transient 
behavior would require models of drastically greater complexity which are not expressible as 
lumped-parameter elements, whose behavior may or may not be expressible in closed form 
analytical functions. In effect, the analysis/simulation problem is an enigmatic trade-off 
between elemental model simplicity and algorithm complexity: a numerically based analysis 
algorithm begs to be discredited and linear while thermal systems are unquestionably contin- 
uous and inevitably non-linear. If element models and their interactive behavior could be 
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linearized and discredited, then all processes could be simply reduced to differences amongst 
linear functions of state variables. However, it is well known from the study of thermody- 
namics that most real processes are path-dependent, rendering such an approach ineffective. 

Future study may investigate similar problems in different fields. For example, the 
finite element analysis of plastic flow and failure of a solid is a problem in continuous, non- 
linear behavior that has been discredited. However, it differs from the thermal power prob- 
lem in that a power system may have a various type of constitutive components while a steel 
truss is constitutively homogeneous throughout its system. 

Another approach, as a corollary to the CV method, may be to impose equilibrium at 
the elements and interpret in some way the mismatch at connecting nodes (the CV approach 
imposed equilibrium at the nodes and allowed the elements to achieve a local equilibrium). It 
has been suggested that an initial-value approach be investigated.. .the number of potential 
paths is still unknown. 

This author wishes to illustrate that this is one approach to the problem, an approach 
which sought a delicate balance between element simplification and algorithmic complexity. 
At the expense of sounding reflective, the author feels compelled to share these uncertainties, 
questions and apprehensions for the edification of future investigators to keep this line of 
research viable. 
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Appendix 1 Recursion Relations for Generating the {R} Vector 

As in the assembly of [K‘], two topology matrices are used to generate { R } : thermal 
topology to generate 1st Law and temperature constitutive terms, and mass topology to gen- 
erate the pressure characteristics and mass continuity. 

1st Law: 



The upper summation limit "n" corresponds to the number of nodes in the system. 
Thermal characteristic constitutive relation: 

For the expander: 



n 



Ri = I m J h J Top T [i,J] 



n 



Ri = I {hj -r\,h JS }7op T [i,j] 

7 = 1 



- fl\ 

R 2 = y Ah -tiTi TTodJ/./I - [Top T [i,j] + 1} 





For the compressor, the usual substitutions apply. 



Pressure characteristic constitutive relation: 



For compressors: 
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(expander) 



«<= 

j =■ 1*- J & -I 



For fluid flow with momentum loss (heat exchangers): 



R,= I 

7 = 1 



p r 



32/7 

v ttD 3 y 






Top m [/,/] 



Continuity constitutive relation: 



*i= im,Top„[i,j] 

7 = 1 
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Appendix 2 Terms of the [K l ] matrix 

The algebraic system described in Chapter 4 must be manipulated to conform to the 
numerical requirements of the [K‘]{Ax}={R} solution methodology. 

Partial derivatives of 1st Law Relations 

As mentioned in Section 5.2, 1st law relations that have environmental interactions are 
not incorporated in K‘. Thus, a first law relation for an adiabatic device with no work transfer 
can be written 



R i (x J )=lm J h J = 0 

j = i 

where p represents the number of fluid streams entering or leaving the control volume for 
device "i". P, v and m at each interconnection are the independent global variables. Hence, 
for each residual relation, the following partial derivatives are of interest: 

dR t dR, dR, 

dP/ dVj’ dm j 



For the 1st Law residual, the previous equations become: 



dR, p . dhj 
■r— -= L m:z— 

dP, j = i J dPj 



dRi 

d^ 



I 



j = i 




dRi p 
TT-= X h 
dm j j = \ J 
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These recursive relation will be combined with the topology matrices in Chapter 7 to yield 
the actual K l entries. 

Partial Derivatives of Characteristic Relations 
Expansion engine characteristic functions: 



(h l -h 2 )-r\ l (h u -h 2s ) = 0 
h x -h u = 0 

s 2t -s l = 0 

Direct differentiation with respect to the independent variables leads to the linearized charac- 
teristic functions: 



dR x dh x 



dh 



dh 2 , ds 



dF, ~ dF x vi ~ ^ dF x v ' ^ ds* Pl 3 p , V1 



dR x dh 2 

dF 2 = ~dF 2 ]v2 



dR x dh x dh x dh * 3 % 

0V! " a^ 1 f 1 " T1, 3v, 1 p 1 + T1, 1 ] p 1 



dR x dh 2 

d7 2 = ~d^ 2 p2 



a /? 2 a/i, dh u 

— -= — -] — -l 

ap, ap, vi ap , 1 
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1 



dR 2 dh x 
dv, dv J p 1 



dh 






dv. 



pi 



dR 2 a/? 2 
dP 2 9v 2 



dRi dSi, ds, 

wr wytf- 



a/? 3 _a% | _a*. 

9v 1 dv, />1 dv, />1 



a /? 3 _ a /? 3 

dP 2 dv 2 



The linearized constitutive relations for the compressor are similar to those for the expander 
engine and an additional characteristic equation for the pressure elevation term: 



R 4 — P, r E P 2 — 0 

dR 4 dR 4 

dP { dP 2 E 



Regenerative heat exchanger characteristic equations: 



Rl=Pl-P 2 - 



16/7 

TtD 3 



m i(v, + v 2 ) = 0 



12 
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r 2 = p 3 -p 4 - 



16/7 

nD' : 



rii](v 3 + v,) = 0 



/34 



/?3 = m l h l -rh 2 h 2 - U R A R In' 



(r 1 -r 4 )-(r 2 -r 3 ) [ 
(r 1 -r 4 )/(r 2 -r 3 ) J 



The partial derivatives are: 



377, _ dR x 

dPT~ l dP 2 



dR x 

drill 




\ 

rii x {v x +v 2 ) 

J 



dR x 

3v, 



16/7^1 2 

r m 2 ( 1 +v 2 ) 

V 7l D 3 / 1 2 



a/?, 

3 v 2 




\ 

m 2 (v, + l) 
J 



dR 2 dR 2 

dP 4 = ~ l dP 3 



dR 2 
5m 3 




\ 

m 3 (v 3 + v 4 ) 

J 



dR 2 
dv 4 




\ 

m 3 (l+v 3 ) 

J 
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dR 2 
9v 3 



^ 1 6/7 

v 7lZ) 3 



m 3 (l +v 4 ) 

J 



Special attention must be paid to the linearization of the rate equation. Obviously, 



dR 3 dR 3 

drhi 9m 2 



Partial derivatives with respect to P and v yield, complex functions that may be replaced by a 
numerical derivative. To implement this, the chain rule is first applied: 

R 3 = g(h i (P,v),T i (P,v),rii i ) 

dR 3 _dg_dh L dg dT) 
dPi ~ dhi dPi v ‘ + 97 ; dPi v ‘ 



dR 3 _dg_dh^ dg dT. 
d Vi ~dh i dv i ir ‘ + dT i dv i lF ‘ 



dg_ 

dh{ 



-ni: 



note also that 



dhi 



and dhi 
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are constitutive relations of the working substance and would exist as subroutines external to 
the [K‘] system, evaluated when called for during the assembly operation. This applies as 
well for 



d J± 

dP, 



and 



d l± 

3v, 



The function dg/dT t shall be evaluated numerically to avoid the analytical function involving 

the logarithm. In general, the numerical derivative for a function of multiple variables with 
respect to the variable x ; is 

dg (xi,x 2 , . . „x„ . . .,x„) g(x i,x 2> . . „x, + Ax,-, . ..,x n ) - g (x„x 2 , . . .,x,-, . . 
dx i Ax, 



For this numerical representation, the contribution to the linearization function from the m h i 

terms subtract to zero since h is not an explicit function of T in this postulation. This leaves 
only the non-linear (logarithmic) portion of the rate equation to be evaluated, denoted f(Tj): 



f(T l )=U AH 



(T, -T 4 )-(T 2 -T 3 ) | 

(r,-r 4 )/(r 2 -r 3 ) J 



The partial derivative of the rate equation residual with respect to (say) P, is 

d_R 3 _ /(T,) _ /(T 1 (Pi + AP„v 1 ),T 2 ,T 3 ,T 4 )-/(r i (P 1 ,v 1 ),T 2 ,r 3 ,T 4 ) 

dP, dPi A P, 



Other partial derivatives are straight forward. 
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The partial derivative relations for the ambient heat exchangers are similar to the 
results for the regenerative heat exchangers with the expected variable and parameter substi- 
tutions. 

Expansion valve characteristic equation: 

R i = h l (P 1 ,v l )-h 2 (P 2 ,v 2 ) 



The partial derivatives are straight forward: 



dR x _ dhi dR , _ 9/z, 

dP x dP x ]v ' 0V! 



a/?,_ dh 2 a/?,_ dh 2 

dp 2 ~~dp} v * d7 2 ~~d7 2 p > 



Mass continuity relations. 



m i -m i + l = 0 



Since m, are independent variables, 



dP,~dv~° 

dR j +1 (in-flow) ) 
drhi l-l(out-flow)J 
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Appendix 3 Recursive Relations for Generation of the [K‘] Matrix 

The partial derivatives of Section 6.2 are combined with the topology matrices to form 
the inner matrices, which are then used to generate the global [K‘] matrix. Each row within 
an inner matrix corresponds to a partial derivative of a constitutive relation for a storage ele- 
ment. 



The recursives for the 1st law interactions are : 



J dPj 



Top T [i,j] 





Top r [/',/| 



(K'jh = hjT°p T [i ,j] 



The recursive relations for the characteristic equations are not uniform between the different 
devices. For the work transfer devices. 



, t 3/t, dh is 



+ Topj.fi,>] 



ds j + 



Topj.fi ,>] 



Topj.fi, j] 



j+Topj.[i,f] 



dP, 



U (Top r [/,/| + 1}} 



Referring to Section 4.3, the recursion for the compressor engine is logically extended by 
replacing l) with (1 — rj c ) . 

For the expansion valve, the temperature characteristic recursion is 
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Top r [/,;'] 



(K‘/ P = dh ; 



dP, 



(K‘j) T v = dhj 



dVj 



Top r [/,y] 



<*;>.« o 



For the regenerative heat exchanger, the partial derivative of the temperature characteristics 
become 



(K‘ j )l = h J To VT [i,j] 



(K'X = df(T t ) 

dP, 



Top r [/,y] 



df(T t ) 

where is the numerical derivative defined in Appendix 1. 



-T7- exists as a subroutine 

BPj 



with standard arguments (T(Pj, Vj)). When the topological argument is non-zero, the corre- 
sponding nodal variable will be passed to the subroutine computing the derivative and appro- 
priately incremented (APj, or AVj) in both the logarithmic function in the numerator and in the 
denominator. Similarly, 






3v, 



Top r [/,y] 
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Here again, the thermal topology is used. 



Pressure characteristic recursions are generated by combining the partial derivatives in 
Appendix 1 with the mass topology. The pressure relations for work transfer devices require 
some modification to the constitutive relation to make it compatible as a numerically gener- 
ated partial derivative. For pressure elevating devices, the pressure elevation is always speci- 
fied on the outlet term, which indicates that indicates that K'y must be able to distinguish inlet 
from outlet. This information lies in the topology matrix where all exiting flows have 
arithmatic value (-1). To implement this, the pressure residual must be re-arranged as fol- 
lows: 



^?i — P x r E P 2 — 0 




— (Vfe) P i — (V^e) P 2 




— 5 To Pm[‘ ;Top_[i ,j] 

r * 
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< 



and 



Top J 1,2 ]=(K l J f , 



^l_j -jTopji.a 

The recursive relation for a pressure characteristic can be generalized to 

< if s^ = |r7 T0P '"'' , } To P-[‘''.'l 

Application to compression devices is simply extended by replacing r E with r c . The recur- 
sions for pressure characteristics experiencing momentum losses are 

(K‘/ p = TopJ/j] 



-32fl \ . 









^ -16/7 ^ 

v 7lD 3 y 



™](i + V T o P j/,,i)fTopj/,y] 



The mass topology matrix is used to modify the variable index on the specific volume so that 
the recursion will generate the exact function. Section 6.3 showed that when 



r=p 1 -p 2 



V 7lD 3 y 



(Wi(v, + v 2 ) 



dR_ 

dv 2 



yirt 3 j 



m\{\ +Vj) 
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; 



The variable index "2", which is the linearizing variable, disappears, but the variable index 
"1" remains. The matrix variable index is j=2, so j=2 must be modified to ensure that the 
variable index "1" appears in the function. This is done by using the arithmatic value of the 
mass topology matrix to modify the variable index since v 2 (an outlet variable) has arithmatic 
value (-1). 

Finally, the mass continuity interaction region can be developed. Clearly, 



(^);=(^)>o 



The mass flux partial derivative for the continuity interaction are identically the mass topol- 
ogy matrix: 



(/g; = Top m [/,y] 



where "i" corresponds to a mass continuity constitutive relation. 
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Appendix 4 Appendix 4: Numerical Solution of [K‘(Xi)]{Ax}={R(Xj)} 

Numerical solution (or reduction) of the [K'(x i k )](Ax k } = (R(x i lt ))=0 system is accom- 
plished using a Newton-Raphson scheme. The index "k" represents the "kth" numerical value 
based on thermodynamic state (x k ). The reduction operation consists of an initial estimate to 
the independent variable vector {x 0 }. (x°) is then used to assign a numerical value to 
[K‘(x 0 )]). The product [K t,0 ]{Ax°) is formed which is the first estimate for the system’s con- 
stitutive relations. The same value of (x°) is used to compute the residual relations (R(x°)}, 
which is then compared to {e}, the error vector. Typically, {R(x°) }>{£}, which requires 
modifying (iterating) the vector (x) until (R(x°) }<{£}. 

The {x} vector is initialized to {x 0 } by the user based on some educated or experience- 
based guess. When {x°) is established, [K‘(x)] assumes its first numerical approximation, 
[K 1,0 ]. Each time the (R(x k )} exceeds the error vector { £ ) , (x) can be updated by 

{x k+l }=-[K'(x k )]-' {*(**)} +{x k } 



{x k + '}={Ax k }+{x k } 



This is simply Newton’s Method solved for (x k+1 ), the improved estimate. 

The error vector {£} is used to check for convergence toward the equilibrium state 
(x^} by comparison with the (R(x k )}. Therefore, {£} should reflect the scale of the constitu- 
tive relation to which it refers. If ( £ ) is suitably chosen, it will ensure quick convergence 
with acceptable accuracy; (£} too large will yield rapid convergence but poor accuracy while 
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on the other hand { e } too small will lead to longer convergence times while generating accu- 
racy that may be excessive for the needs of the analysis. The method proven useful in Dr. 
Russo’s CAT methodology is to use a fractional error based on the constitutive relation of 
interest: 



[R,(x°)[ 

e ' IT" 

where "n" represents the number of independent variables. For a small system, n may be too 
small to render desired accuracy, in which case some other suitable divisor may be used. This 
is a proven, simple and consistent system scaling technique. 

When {x^} is attained, the post processing calculations are executed since all the inde- 
pendent variables are known and all the post-processing relations are functions of the inde- 
pendent variables. As a minimum, the 1st Law relations for non-zero heat and work transfer 
are computed to complete the system analysis. Additional computations based on the now 
determined state variables may now be conducted as specified by the designer/operator. 



72 




Thesis 

£6766 Stanton 

C *1 An algorithm for con- 

trol volume analysis of 
cryogenic systems. 



COtl- 
s of 



